
clc; clear; close all;

%% 1.初始化
N = 20000;           % 总迭代步数
Transient = 1000;    % 扔掉的初始过渡步数(消除瞬态)

% 论文中常用的超混沌参数
a = 0.8;   % QO-IM
b = 1.5;   % SO-IM
c = 1.7;   % TO-GM
d = 1.8;   % CO-GM

% 论文中的初始状态
x1_0 = 0.5;  phi1_0 = 0.5;   % QO-IM
x2_0 = 0.1;  phi2_0 = 1.0;   % SO-IM
x3_0 = 0.5;  phi3_0 = -1.0;  % TO-GM
x4_0 = 0.5;  phi4_0 = 1.0;   % CO-GM

%% 2. 迭代计算
% 分别存储四个映射的 (x, phi)
[x_QO, phi_QO] = QO_IM_map(a, x1_0, phi1_0, N, Transient);
[x_SO, phi_SO] = SO_IM_map(b, x2_0, phi2_0, N, Transient);
[x_TO, phi_TO] = TO_GM_map(c, x3_0, phi3_0, N, Transient);
[x_CO, phi_CO] = CO_GM_map(d, x4_0, phi4_0, N, Transient);

%% 3. 作图
figure('Color','w');

% (a) QO-IM
subplot(2,2,1);
plot(x_QO, phi_QO, '.', 'Color',[0.8,0.6,0.1], 'MarkerSize',2);
xlabel('\itx_1'); ylabel('\phi_1');
title('(a) QO-IM');
axis tight; box on; grid on;

% (b) SO-IM
subplot(2,2,2);
plot(x_SO, phi_SO, '.', 'Color',[0,0.6,0.4], 'MarkerSize',2);
xlabel('\itx_2'); ylabel('\phi_2');
title('(b) SO-IM');
axis tight; box on; grid on;

% (c) TO-GM
subplot(2,2,3);
plot(x_TO, phi_TO, '.', 'Color',[0.9,0.3,0.2], 'MarkerSize',2);
xlabel('\itx_3'); ylabel('\phi_3');
title('(c) TO-GM');
axis tight; box on; grid on;

% (d) CO-GM
subplot(2,2,4);
plot(x_CO, phi_CO, '.', 'Color',[1.0,0.5,0.0], 'MarkerSize',2);
xlabel('\itx_4'); ylabel('\phi_4');
title('(d) CO-GM');
axis tight; box on; grid on;

sgtitle('四种离散忆阻振荡超混沌映射的吸引子');

%% 超混沌映射的封装函数

%% 1. QO-IM
%    x1(n+1) = 0.2(-1)^n + a( (phi1(n))^2 - 2 ) * x1(n)
%    phi1(n+1) = phi1(n) + x1(n)
function [x_data, phi_data] = QO_IM_map(a, x0, phi0, N, Transient)
    x = x0;    phi = phi0;
    x_data = zeros(1,N-Transient);
    phi_data = zeros(1,N-Transient);

    for n = 1:N
        x_next = 0.2*((-1)^n) + a*((phi^2)-2)*x;
        phi_next = phi + x;
        
        x = x_next;
        phi = phi_next;

        if n > Transient
            idx = n - Transient;
            x_data(idx) = x;
            phi_data(idx) = phi;
        end
    end
end

%% 2. SO-IM
%    x2(n+1) = 0.2(-1)^n + b sin(pi phi2(n)) * x2(n)
%    phi2(n+1) = phi2(n) + x2(n)
function [x_data, phi_data] = SO_IM_map(b, x0, phi0, N, Transient)
    x = x0;   phi = phi0;
    x_data = zeros(1,N-Transient);
    phi_data = zeros(1,N-Transient);

    for n = 1:N
        x_next = 0.2*((-1)^n) + b*sin(pi*phi)*x;
        phi_next = phi + x;

        x = x_next;
        phi = phi_next;

        if n > Transient
            idx = n - Transient;
            x_data(idx) = x;
            phi_data(idx) = phi;
        end
    end
end

%% 3. TO-GM
%    x3(n+1) = (-1)^n + c sin(phi3(n)) * x3(n)
%    phi3(n+1)=0.86phi3(n)+0.12phi3(n)|phi3(n)|-0.02phi3(n)^3+0.1 x3(n)
function [x_data, phi_data] = TO_GM_map(c, x0, phi0, N, Transient)
    x = x0;   phi = phi0;
    x_data = zeros(1,N-Transient);
    phi_data = zeros(1,N-Transient);

    for n = 1:N
        x_next = ((-1)^n) + c*sin(phi)*x;
        phi_next = 0.86*phi + 0.12*phi*abs(phi) - 0.02*(phi^3) + 0.1*x;

        x = x_next;
        phi = phi_next;

        if n > Transient
            idx = n - Transient;
            x_data(idx) = x;
            phi_data(idx) = phi;
        end
    end
end

%% 4. CO-GM
%    x4(n+1) = 0.8(-1)^n + d cos(phi4(n)) * x4(n)
%    phi4(n+1)=0.1 tanh(phi4(n)^3) + 0.99 phi4(n) + 0.1 x4(n)
function [x_data, phi_data] = CO_GM_map(d, x0, phi0, N, Transient)
    x = x0;   phi = phi0;
    x_data = zeros(1,N-Transient);
    phi_data = zeros(1,N-Transient);

    for n = 1:N
        x_next = 0.8*((-1)^n) + d*cos(phi)*x;
        phi_next = 0.1*tanh(phi^3) + 0.99*phi + 0.1*x;

        x = x_next;
        phi = phi_next;

        if n > Transient
            idx = n - Transient;
            x_data(idx) = x;
            phi_data(idx) = phi;
        end
    end
end
